#Install packages if needed
remotes::install_github("eliocamp/ggnewscale@dev")
remotes::install_github("YuLab-SMU/ggtree")
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("Biostrings")
BiocManager::install("ggtree")
#For document rendering
tinytex::install_tinytex()
# general
library(tidyverse)
library(readxl)
library(plyr)
library(dplyr)
library(here)
"%notin%" = Negate("%in%")
here::here()
library(tinytex)
# graphics
library(ggplot2)
library(ggnewscale)
library(viridis)
# phylogenetic & other tools
library(ape)
library(Biostrings)
library(ggtree)
library(rotl) #for phylogenetic analyses, get all the species? from Hinchliff et al. 2015 PNAS
library(stringr)
library(tidytree)
NOTES:
Using combination of my trial code and Matt Savoca and company’s code to build a phylo tree with associated trait and %FO data.
A lot of troubleshooting later and I found that the species list needs to go on the far left side of the data, all other data/factors need to be added to the RHS of the species list
# data contains sp list, class, order & family + traits.
# Probable life stage data
my_prey_prob = read.csv(here("./data/output_data/prey_trait_select.csv"), header = TRUE) %>%
dplyr::select(prey_sp, prey_class:prey_family, life_stage:maxM, -c(diet_likely,
life_note)) %>%
filter(prey_sp %notin% c("Lampanyctus mexicanus", "Janthina exigua"))
# Need to remove these two species as they cause issues with ROTL due to naming
# 306 species and 13 variables of interest excluding prey_sp str(my_prey_prob)
Most of the variables except the diet metrics need to be factors.
sapply(my_prey_prob, class)
## prey_sp prey_class prey_order prey_family
## "character" "character" "character" "character"
## life_stage vert_habitat horz_habitat diel_migrant
## "character" "character" "character" "numeric"
## diel_migrant_cat refuge_cat season_migrant season_cat
## "character" "character" "numeric" "character"
## gregarious_primary maxFO maxN maxM
## "character" "numeric" "numeric" "numeric"
my_prey_prob[, 2:13] <- lapply(my_prey_prob[, 2:13], factor)
my_prey_prob$maxFO[my_prey_prob$maxFO == 0] <- NA
my_prey_prob$maxN[my_prey_prob$maxN == 0] <- NA
my_prey_prob$maxM[my_prey_prob$maxM == 0] <- NA
# Tree build 1
breaks <- c(seq(1, nrow(my_prey_prob), 50), nrow(my_prey_prob) + 1) # why are we doing this?
# looking up each of the species in dataset to a phylogeny, doing it by sets of
# 50 if there's an NA (species that didn't match a value in known phylogeny),
# breaks the whole chunk of 50 species
for (i in 1:(length(breaks) - 1)) {
taxa <- as.character(my_prey_prob$prey_sp[breaks[i]:(breaks[i + 1] - 1)])
taxa <- taxa[taxa != "" & !is.na(taxa)]
resolved_namest <- tnrs_match_names(taxa) # I think this is where all the extra species fall out
resolved_namest <- resolved_namest[!is.na(resolved_namest$unique_name), ] # ignore an NA
if (i == 1) {
resolved_namess <- resolved_namest
} else {
resolved_namess <- rbind(resolved_namess, resolved_namest)
}
}
resolved_names <- resolved_namess
resolved_names <- resolved_names[resolved_names$flags != "INCERTAE_SEDIS_INHERITED",
]
# original tree based on simple taxon datset
my_tree_prob <- tol_induced_subtree(ott_ids = resolved_names$ott_id, label_format = "name")
my_tree_prob$tip.label <- gsub("_", " ", my_tree_prob$tip.label) # removes underscore between genus and species names
my_tree_prob$tip.label <- str_extract(my_tree_prob$tip.label, "[A-Z][a-z]+ [a-z]+")
my_tree_prob <- compute.brlen(my_tree_prob, method = "Grafen", power = 1/2) #add branch lengths to my tree using the Grafen (1989) method
my_tree_prob <- ladderize(my_tree_prob, right = TRUE)
View(my_tree_prob)
write.tree(my_tree_prob, here("./data/output_data/albacore_diet_tree_prob"))
my_tree_prob <- read.tree(here("./data/output_data/albacore_diet_tree_prob"))
#306 species
NOTE: Here we edit the tree tip labels - species names - for my_tree_prob and check data dimensions line up with original df containing trait info my_prey_prob.
# Edit tip labels
my_tree_prob$tip.label <- as.factor(sub("_", " ", my_tree_prob$tip.label))
# Check that the number of rows matches the number of tip labels
nrow(my_prey_prob) == length(my_tree_prob$tip.label) #TRUE
## [1] TRUE
# Can check this manually also nrow(my_prey_prob) # 306
# length(my_tree_prob$tip.label) # 306
NOTE: We have in the past encountered a data frame dimension issue when plotting our data values onto the tree. There are also several species names that do not line up. This issue is addressed in the chunk below.
# my_prey_prob$prey_sp <- gsub(' ', '_', my_prey_prob$prey_sp) #for ease of
# plotting take away ' '
tree_names_prob = data.frame(sort(my_tree_prob$tip.label)) #sort the species names from the phylo tree
# There are 306 sort brings tree_names_prob from 301 to 300 because of an NA in
# the tip labels
prey_names_prob = data.frame(sort(my_prey_prob$prey_sp)) #sort the species names from the original list
nrow(tree_names_prob)
## [1] 306
nrow(prey_names_prob) #Needs to be 306 -- 306
## [1] 306
# Sometimes one species didn't make it to the tree
NOTE: Here we identify which species names do not match from our original species and trait df, compared to that which parsed through phylogenetic info from ROTL project. These were:
Original data my_prey_prob, and checked via prey_names_prob: Diplodus sargus, Leuroglossus stilbius, Liocranchia reinhardtii, Mullus barbatus barbatus, Notoscopelus kroyeri, Robustosergia robusta.
Tree data my_tree_prob, and checked via tree_names_prob: Bathylagus stilbius, Liocranchia reinhardti, Mullus barbatus, NA, Notoscopelus elongatus, Sergia robusta.
tree_names_prob == prey_names_prob
## here we check the names that are in the prey data but not on the tip labels
## of the tree
names_not_in_tree_prob = prey_names_prob %>%
filter(sort.my_prey_prob.prey_sp. %notin% tree_names_prob$sort.my_tree_prob.tip.label.)
print(names_not_in_tree_prob$sort.my_prey_prob.prey_sp.)
## [1] "Diplodus sargus" "Leuroglossus stilbius"
## [3] "Liocranchia reinhardtii" "Mullus barbatus barbatus"
## [5] "Notoscopelus kroyeri" "Robustosergia robusta"
## check names that are in the tree tip labels but not in the prey data
names_not_in_prey_prob = tree_names_prob %>%
filter(sort.my_tree_prob.tip.label. %notin% prey_names_prob$sort.my_prey_prob.prey_sp.)
print(names_not_in_prey_prob$sort.my_tree_prob.tip.label.)
## [1] Bathylagus stilbius Liocranchia reinhardti Mullus barbatus
## [4] NA Notoscopelus elongatus Sergia robusta
## 306 Levels: Abralia redfieldi Abraliopsis affinis ... Walvisteuthis rancureli
NOTE: These taxa need to be re-aligned via detailed code below.
## here we decide to match the names in the prey data to the tip labels, so we create a dataset with the names in the prey data that are in the tree (I see this looks kind of confusing with the double negative, the names not not in the tree). These are the names that do not need to be fixed, which we will bind the fixed names to afterwards so we do not get doubles of the same names (with different spelling).
my_prey_keep_prob = my_prey_prob %>%
filter(prey_sp %notin% names_not_in_tree_prob$sort.my_prey_prob.prey_sp.)
#300 names to keep, and 6 need to be corrected
## these are the names we are going to fix to match the tip labels of the tree, given they are synonyms/misspellings. These are the same names in 'names_not_in_tree_adult', except they now have the columns of data which will allow use to merge it to 'my_prey_keep_adult' after we fix the names.
my_prey_fix_prob = my_prey_prob %>%
filter(prey_sp %in% names_not_in_tree_prob$sort.my_prey_prob.prey_sp.)
#6 names to fix
## this step isn't necessary, but it does make it easier in the names are in alphabetical order
my_prey_fix_prob = my_prey_fix_prob[order(my_prey_fix_prob$prey_sp),]
# the names have to be characters for the name reassignment to work, it will give an error if we don't do this (they are factors before this)
my_prey_fix_prob$prey_sp <- as.character(my_prey_fix_prob$prey_sp)
## we make these into x2 and y2 because it makes it much easier to write the next section
x2 = my_prey_fix_prob
y2 = as.vector(names_not_in_prey_prob$sort.my_tree_prob.tip.label.)
unique(x2$prey_sp)
## [1] "Diplodus sargus" "Leuroglossus stilbius"
## [3] "Liocranchia reinhardtii" "Mullus barbatus barbatus"
## [5] "Notoscopelus kroyeri" "Robustosergia robusta"
y2
## [1] "Bathylagus stilbius" "Liocranchia reinhardti" "Mullus barbatus"
## [4] "NA" "Notoscopelus elongatus" "Sergia robusta"
#here we look at the names in both of the lists and see if there are misspellings/synonyms in the names that caused the discrepancy between the name lists. As you can see this is true for all of the species in these lists expect one. 'Diplodus sargus' does not have a match in the tree, since the only other tip labels left in the tree is an NA after the rest of the reassignments.
# if you want to see how this works you can run individual pieces of this before running the whole thing. For example x2[4,1] is "Leuroglossus stilbius", and y2[1] is "Bathylagus stilbius". Then after running this code we match the name in the prey data (x2) to the names in the tree (y2), so "Leuroglossus stilbius" becomes "Bathylagus stilbius".
x2[1,1] = y2[4]; #--> Diplodus sargus = NA
x2[2,1] = y2[1]; #--> Leuroglossus stilbius = Bathylagus stilbius
x2[3,1] = y2[2]; #--> Liocranchia reinhardtii = Liocranchia reinhardti
x2[4,1] = y2[3]; #--> Mullus barbatus barbatus = Mullus barbatus
x2[5,1] = y2[5]; #--> Notoscopelus kroyeri = Notoscopelus elongatus
x2[6,1] = y2[6]; #--> Robustosergia robusta = Sergia robusta
# we need to do the opposite for one of the names here "Diplodus sargus", which is in the prey data but not in the tip labels. So, we are replacing the NA in the tip labels
tip_labels_prob <- as.vector(my_tree_prob$tip.label)
tip_labels_prob[28] ##this is the NA we want to replace
## [1] "NA"
tip_labels_prob[28] = "Diplodus sargus"
## replacing the tip labels with the additional name
my_tree_prob$tip.label <- tip_labels_prob #Still have 306 tip labels
# here we bind the fixed names to the names that didnt need to be fixed
my_prey_prob = rbind(my_prey_keep_prob, x2)
#here we get rid of any of the names in the data that didn't have a match in the tree tip labels, which was only 'Diplodus sargus'
#--> This is where one species drops out!
# Do not run this bit
#my_prey_prob = my_prey_prob %>%
#--> This is where one species drops out!
# filter(prey_sp %in% my_tree_prob$tip.label)
nrow(my_prey_prob);length(my_tree_prob$tip.label) ## one of those tip labels is an NA
## [1] 306
## [1] 306
#the prey names need to be the rownames for the heatmap to work on the tree
#rownames(my_prey_prob) = NULL #this did not solve the row names thing
rownames(my_prey_prob) = my_prey_prob$prey_sp
#tip labels need to be characters to pipe onto the tree
my_tree_prob$tip.label <- as.character(my_tree_prob$tip.label)
str(my_tree_prob)
## List of 5
## $ edge : int [1:578, 1:2] 307 308 309 310 311 312 313 314 315 316 ...
## $ edge.length: num [1:578] 0.00164 0.24959 0.0066 0.01789 0.03479 ...
## $ Nnode : int 273
## $ node.label : chr [1:273] "mrcaott42ott150" "mrcaott42ott49" "mrcaott42ott658" "Clupeocephala" ...
## $ tip.label : chr [1:306] "Sebastes wilsoni" "Sebastes zacentrus" "Sebastes proriger" "Sebastes brevispinis" ...
## - attr(*, "class")= chr "phylo"
## - attr(*, "order")= chr "cladewise"
Basic prey species phylogeny graph
# Phylotree without labels
prob_basic = ggtree(my_tree_prob, layout = "circular")
prob_basic
# Phylotree with species labels
prob_basic_lab = ggtree(my_tree_prob, layout = "circular") + geom_tiplab(size = 2)
prob_basic_lab
Class-based phylogeny graph
#Use basic + prey_sp labels + prey_class
# creating a dataset that splits the species data by class
prey_class_prob_info <- split(x = my_prey_prob$prey_sp, f = my_prey_prob$prey_class)
unique(my_prey_prob$prey_class)
## [1] Actinopterygii Cephalopoda Gastropoda Hexanauplia Hydrozoa
## [6] Malacostraca Thaliacea
## 7 Levels: Actinopterygii Cephalopoda Gastropoda Hexanauplia ... Thaliacea
# splitting the tree species by class
my_tree_prob_class <- groupOTU(my_tree_prob, prey_class_prob_info)
as.character(sort(unique(my_prey_prob$prey_class)))
## [1] "Actinopterygii" "Cephalopoda" "Gastropoda" "Hexanauplia"
## [5] "Hydrozoa" "Malacostraca" "Thaliacea"
#creating a palette for the class split
#pal_prob_class = c("grey40", '#0047ab', '#751308', '#4B0082', '#F05E23', '#013220', '#B80F0A', "#66C2A5", "#3288BD") ## the grey40 is needed for a branch between branches, but doesn't assign to a class
pal_prob_class = c('#0047ab', #Actinopterygii
#'#751308', #Branchiopoda #older teal colour "#66C2A5" #Not in prob data just was used for adult data
'#4B0082', #Cephalopoda
"#AD15E2", #Gastropoda #old orange '#F05E23'
'#B80F0A', #Hexanauplia
"#773405", #Hydrozoa #'#013220'
"#E86103", #Malacostraca '
"#3288BD", #Thaliacea
"grey40" ## the grey40 is needed for a branch between branches, but doesn't assign to a class
)
## plotting the tree without tip labels
prob_class <- ggtree(my_tree_prob_class, aes(color = group), layout = 'circular') +
theme(legend.position = c(0.5,0.50),
legend.title = element_text(size = 18),
legend.text = element_text(size = 18))+
scale_colour_manual('Class', aesthetics = c('colour'), values = pal_prob_class,
breaks = c("Actinopterygii", "Cephalopoda","Gastropoda", "Hexanauplia", "Hydrozoa", "Malacostraca","Thaliacea"),
labels = c("Actinopterygii", "Cephalopoda","Gastropoda", "Hexanauplia", "Hydrozoa", "Malacostraca","Thaliacea"))
prob_class
# plotting the tree with tip labels
prob_class_lab <- ggtree(my_tree_prob_class, aes(color = group), layout = 'circular') +
theme(legend.position = c(0.5,0.50),
legend.title = element_text(size = 18),
legend.text = element_text(size = 18))+
scale_colour_manual('Class', aesthetics = c('colour', 'fill'), values = pal_prob_class,
breaks = c("Actinopterygii", "Cephalopoda","Gastropoda", "Hexanauplia", "Hydrozoa", "Malacostraca","Thaliacea"),
labels = c("Actinopterygii", "Cephalopoda","Gastropoda", "Hexanauplia", "Hydrozoa", "Malacostraca","Thaliacea")) +
geom_tiplab(size = 2)
prob_class_lab
Frequency Occurrence
#Overlay maxFO (frequency of occurrence data) on a basic prey tree coloured by Class
#Use prey_class, maxFO
prob_maxfo <- gheatmap(prob_class, my_prey_prob[,"maxFO", drop = FALSE],
offset = 0, width = 0.05, colnames = FALSE) +
theme(#legend.position = c(0.5,0.50),
legend.title = element_text(size = 18),
legend.text = element_text(size = 18))+
scale_fill_viridis_c(name = "Percent frequency (%FO)", na.value = 'grey')
prob_maxfo
Abundance
#Overlay maxN (abundance data) on a basic prey tree coloured by Class
#Use prey_class, maxN
my_prey_prob$maxN <- as.numeric(my_prey_prob$maxN)
prob_maxn <- gheatmap(prob_class, my_prey_prob[,"maxN", drop = FALSE],
offset = 0, width = 0.05, colnames = FALSE) +
theme(#legend.position = c(0.5,0.50),
legend.title = element_text(size = 18),
legend.text = element_text(size = 18))+
scale_fill_viridis_c(name = "Percent abundance (%N)", na.value = 'grey')
prob_maxn
Biomass
#Overlay maxM (biomass data) on a basic prey tree coloured by Class
#Use prey_class, maxM
prob_maxm <- gheatmap(prob_class, my_prey_prob[,"maxM", drop = FALSE],
offset = 0, width = 0.05, colnames = FALSE) +
theme(#legend.position = c(0.5,0.50),
legend.title = element_text(size = 18),
legend.text = element_text(size = 18))+
scale_fill_viridis_c(name = "Percent biomass (%KG)", na.value = 'grey')
prob_maxm
Habitat Use Traits
For habitat use traits: we are going to plot multiple rings of data around the phylogenetic trees using vert_habitat, horz_habitat, diel_migrant (change to yes/no), and season_migrant (yes/no)
#creating this plot one iteration at a time. Adding one layer of the heatmap, then allowing for another legend on the plot, and then another layer of the heatmap
habitat_traits_prob1 <- gheatmap(prob_basic, my_prey_prob[, 'vert_habitat',drop=FALSE],
offset= 0, width=0.05, colnames = FALSE) +
scale_fill_viridis_d(name = "Vertical Habitat (VH)",
direction = -1,
breaks = c("benthic", "demersal", "epipelagic", "mesopelagic", "bathypelagic"),
limits = c("benthic", "demersal", "epipelagic", "mesopelagic", "bathypelagic"),
na.translate = TRUE,
option = 'D',
na.value = 'grey')+
theme(legend.position = c(1,0.80),
legend.title = element_text(size = 18),
legend.text = element_text(size = 18))+
scale_y_continuous(expand = c(0,4))+ #this opens up the gap for the Vertical habitat label around the tree
annotate('text', x = 1.03, y = -7, label = 'VH', angle = -85, size = 4)
habitat_traits_prob1.5 <- habitat_traits_prob1 + new_scale_fill()
habitat_traits_prob2 <- gheatmap(habitat_traits_prob1.5, my_prey_prob[ ,'horz_habitat',drop=FALSE],
offset=0.05, width=0.05, colnames = F) +
scale_fill_viridis_d(name = "Horizontal Habitat (HH)",
option = "D",
direction = -1,
breaks = c("reef-associated", "coastal", "continental shelf","continental slope", "oceanic"),
limits = c("reef-associated", "coastal", "continental shelf","continental slope", "oceanic"),
na.translate = TRUE,
na.value = 'grey')+
theme(legend.position = c(1,0.683),
legend.title = element_text(size = 18),
legend.text = element_text(size = 18))+
scale_y_continuous(expand = c(0,4))+
annotate('text', x = 1.08, y = -7.1, label = 'HH', angle = -85, size = 4)
habitat_traits_prob2.5 <- habitat_traits_prob2 + new_scale_fill()
habitat_traits_prob3 <- gheatmap(habitat_traits_prob2.5, my_prey_prob[ , 'diel_migrant',drop=FALSE], offset=0.10, width=0.05, colnames = F) +
scale_fill_manual(name = "Diel Migrant (DM)",
breaks = c("0", "1"),
limits = c("0", "1"),
labels = c("no","yes"),
values = c("0"="#FDE725FF", "1"="#39568CFF"),
na.value = 'grey')+
theme(legend.position = c(1,0.5985),
legend.title = element_text(size = 18),
legend.text = element_text(size = 18))+
scale_y_continuous(expand = c(0,4))+
annotate('text', x = 1.13, y = -7.5, label = 'DM', angle = -85, size = 4)
habitat_traits_prob3.5 <- habitat_traits_prob3 + new_scale_fill()
habitat_traits_prob_final <- gheatmap(habitat_traits_prob3.5,
my_prey_prob[ ,'season_migrant',drop=FALSE],
offset=0.15, width=0.05, colnames = F) +
theme(legend.position = c(1.05,0.50),
legend.title = element_text(size = 18),
legend.text = element_text(size = 18))+
scale_fill_manual(name = "Seasonal Migrant (SM)",
breaks = c("0", "1"),
limits = c("0", "1"),
labels = c("no", "yes"),
values = c("0"="#FDE725FF", "1"="#39568CFF"),
na.value = 'grey') +
scale_y_continuous(expand = c(0,4))+
annotate('text', x = 1.18, y = -7.5, label = 'SM', angle = -85, size = 4)
#Check graph
habitat_traits_prob_final
Habitat use and aggregation trait
Adding aggregation to habitat use figure
#creating this plot one iteration at a time. Adding one layer of the heatmap, then allowing for another legend on the plot, and then another layer of the heatmap
habitat_traits_prob1 <- gheatmap(prob_basic, my_prey_prob[, 'vert_habitat',drop=FALSE],
offset= 0, width=0.05, colnames = FALSE) +
scale_fill_viridis_d(name = "Vertical Habitat (VH)",
direction = -1,
breaks = c("benthic", "demersal", "epipelagic", "mesopelagic", "bathypelagic"),
limits = c("benthic", "demersal", "epipelagic", "mesopelagic", "bathypelagic"),
na.translate = TRUE,
option = 'D',
na.value = 'grey')+
theme(legend.position = c(1,0.80),
legend.title = element_text(size = 18),
legend.text = element_text(size = 18))+
scale_y_continuous(expand = c(0,4))+ #this opens up the gap for the Vertical habitat label around the tree
annotate('text', x = 1.03, y = -7, label = 'VH', angle = -85, size = 4)
habitat_traits_prob1.5 <- habitat_traits_prob1 + new_scale_fill()
habitat_traits_prob2 <- gheatmap(habitat_traits_prob1.5, my_prey_prob[ ,'horz_habitat',drop=FALSE],
offset=0.05, width=0.05, colnames = F) +
scale_fill_viridis_d(name = "Horizontal Habitat (HH)",
option = "D",
direction = -1,
breaks = c("reef-associated", "coastal", "continental shelf","continental slope", "oceanic"),
limits = c("reef-associated", "coastal", "continental shelf","continental slope", "oceanic"),
na.translate = TRUE,
na.value = 'grey')+
theme(legend.position = c(1,0.683),
legend.title = element_text(size = 18),
legend.text = element_text(size = 18))+
scale_y_continuous(expand = c(0,4))+
annotate('text', x = 1.08, y = -7.1, label = 'HH', angle = -85, size = 4)
habitat_traits_prob2.5 <- habitat_traits_prob2 + new_scale_fill()
habitat_traits_prob3 <- gheatmap(habitat_traits_prob2.5, my_prey_prob[ , 'diel_migrant',drop=FALSE], offset=0.10, width=0.05, colnames = F) +
scale_fill_manual(name = "Diel Migrant (DM)",
breaks = c("0", "1"),
limits = c("0", "1"),
labels = c("no","yes"),
values = c("0"="#FDE725FF", "1"="#39568CFF"),
na.value = 'grey')+
theme(legend.position = c(1,0.5985),
legend.title = element_text(size = 18),
legend.text = element_text(size = 18))+
scale_y_continuous(expand = c(0,4))+
annotate('text', x = 1.13, y = -7.5, label = 'DM', angle = -85, size = 4)
habitat_traits_prob3.5 <- habitat_traits_prob3 + new_scale_fill()
habitat_traits_prob4 <- gheatmap(habitat_traits_prob3.5,
my_prey_prob[ ,'season_migrant',drop=FALSE],
offset=0.15, width=0.05, colnames = F) +
theme(legend.position = c(1.05,0.50),
legend.title = element_text(size = 18),
legend.text = element_text(size = 18))+
scale_fill_manual(name = "Seasonal Migrant (SM)",
breaks = c("0", "1"),
limits = c("0", "1"),
labels = c("no", "yes"),
values = c("0"="#FDE725FF", "1"="#39568CFF"),
na.value = 'grey') +
scale_y_continuous(expand = c(0,4))+
annotate('text', x = 1.18, y = -7.5, label = 'SM', angle = -85, size = 4)
habitat_traits_prob4.5 <- habitat_traits_prob4 + new_scale_fill()
habitat_traits_prob_agg_final <- gheatmap(habitat_traits_prob4.5,
my_prey_prob[ ,"gregarious_primary",drop=FALSE],
offset=0.20, width=0.05, colnames = F) +
theme(legend.position = c(1.05,0.50),
legend.title = element_text(size = 18),
legend.text = element_text(size = 18))+
scale_fill_viridis_d(name = "Gregariousness (GG)",
direction = -1,
breaks = c("solitary", "shoaling","schooling"),
limits = c("solitary", "shoaling","schooling"),
na.translate = TRUE,
option = "D", #'magma',
begin = 0.2,end = 0.8,
guide = guide_legend(order = 1),
na.value = 'grey')+
scale_y_continuous(expand = c(0,4))+
annotate('text', x = 1.23, y = -7.5, label = 'GG', angle = -85, size = 4)
#Check graph
habitat_traits_prob_agg_final
There are numerous issues with saving output either via ggsave() or dev.copy2pdf() functions. I’m batching the output saves below and doing all at once.
Basic phylogeny
Including code for saving this graph. Code for saving output for all other figures below are included in the .Rmd and not in the rendered document.
Phylogeny by Class
Phylogeny by diet use
Phylogeny by habitat Use Traits